!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine names
  !title
  print *,' '
  print *,'              DROPS-FD_MC 1.0'
  print *,' '
  !authors
  print *,'         Dr. Julio C Armas Perez'
  print *,'         Prof. Dr. J.P. Hernandez-Ortiz'
  print *,'         Prof. Dr. O. Guzman'
  print *,'         Prof. Dr. J.J. de Pablo'
  print *,' '
  !afiliation
  print *,'            Department of Materials '
  print *,'     Universidad Nacional de Colombia, Medellin'
  print *,' '
  print *,'             Department of Physics'
  print *,'     Universidad Autonoma Metropolitana, Mexico '
  print *,' '
  print *,'        Institute for Molecular Engineering'
  print *,'             University of Chicago'
  print *,' '
end subroutine names
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
function S(U)
  use rbf_type
  implicit none
  real(dp), intent(in) :: U
  real(dp) S
  S = 0.25d0 + 0.75d0*sqrt( 1.d0 - 8.d0/(3.d0*U))
end function S
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine q_anchoring
  use rbf_type
  use variable, only: normal, nn, chi, nband, width, coverage,tangential,tangential_2,xm
  use variable, only: surfactant, lband, ao_s, Sbulk,planar,ao_w,plan
  implicit none
  integer(i4b) i, j
  real(dp) r, transition
  real(dp), dimension(3):: nrandom, dx
  real(dp), dimension(:,:), allocatable:: qo_s,qo_w,qo_w_2

  allocate( qo_s  (5,nn) ); qo_s = 0.d0
  allocate( qo_w  (5,nn) ); qo_w = 0.d0
  allocate( qo_w_2(5,nn) ); qo_w_2 = 0.d0
  ! Boundary nodes and Prefered anchoring
  do i = 1, nn

     ! UNIFORM CONCENTRATION
     if ( surfactant == 0 )then 
        chi(i) = transition( coverage )
     else if ( surfactant == 1 ) then
        if (xm(3*i).lt.0.d0 )then
           chi(i) = 0.d0 
        else
           chi(i) = 1.d0
        endif
     else if ( surfactant == 2 .or.surfactant == 8 ) then
        if(plan(i))then
           chi(i) = 0.d0
        else
           chi(i) = 1.d0
        endif
     else if ( surfactant == 3 ) then
        chi(i) = 0.d0
   endif


     ! Homeotropic --> Surfactant 
     !  planar  --- > Water
  
     call director_tensor( Sbulk, normal      (3*i - 2:3*i), qo_s  (:,i) )  ! qo_s is the tensor for harmonic term 
     call director_tensor( Sbulk, tangential  (3*i - 2:3*i), qo_w  (:,i) )  ! qo_w is the tensor for harmonic term 
     call director_tensor( Sbulk, tangential_2(3*i - 2:3*i), qo_w_2(:,i) )  ! qo_w is the tensor for harmonic term
  end do
  

  do i = 1, nn
     
     if(surfactant==3)then
        if (xm(3*i).lt.0.d0 )then
           call transformation( qo_w(:,i), ao_w(:,i))
        else
           call transformation( qo_w_2(:,i), ao_w(:,i) )
        endif
     else if(surfactant==2.or.surfactant==8)then
        if(plan(i))then
           call transformation( qo_w(:,i), ao_w(:,i))
        else
           call transformation( qo_s(:,i), ao_s(:,i) )
        endif
     else 
        call transformation( qo_s(:,i), ao_s(:,i) )
        call transformation( qo_w(:,i), ao_w(:,i) )
     endif
  end do
  
  deallocate( qo_s,qo_w )
end subroutine q_anchoring
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine q_initial
  use rbf_type

  use variable, only : Sinitial,nodes,a,xm,direction,nn,nx,ny,nz,nee_cube,pi,Sbulk,planar,normal,chi,q0
  use variable, only : surfactant,lband,tangential,tangential_2,planar,wall,plan,q
  use random_number2
  use ranftest

  implicit none
!  logical(lgt) flag
  integer(i4b) i,j,k,imc
  real(dp) r
  real(dp), dimension(3):: nrandom, dx,dx_norm,dx_plan,dx_tang,dx_tang2,dx_s,dx_w,dx_random
  real(dp), dimension(:,:), allocatable:: qo_s, qo_w
  real(dp) :: thetiux, phiux, omegiux
  real(dp) :: theta, phi, omega,x,y,z,rho,s0
!-----------------------------------------------------------
  
allocate( qo_s(5,nn), qo_w(5,nn) ); qo_s = 0.d0; qo_w = 0.d0

! Boundary nodes and Prefered anchoring

dx = 0.d0; dx(3) = 1.d0

 
do i = 1, nn
     
     ! normal vector
    
     dx_norm  = normal(3*i - 2:3*i)
     dx_tang  = tangential(3*i - 2:3*i) 
     dx_tang2 = tangential_2(3*i - 2:3*i)
     call RANDOM_NUMBER( nrandom )
     nrandom (1) =(2.d0*nrandom (1)-1.d0)
     nrandom (2) =(2.d0*nrandom (2)-1.d0)
     nrandom (3) =0.d0
     r       = sqrt( dot_product( nrandom, nrandom ) )
     dx_random = nrandom/r
    
     if (surfactant.eq.0) then ! the subtrate on the bottom and the top are equal 
        dx_s =  dx_norm
        dx_w =  dx_tang
     else if(surfactant.eq.1)then
        dx_s =  dx_norm
        dx_w =  dx_tang    
     else if(surfactant.eq.2)then   
        dx_s =  dx_norm
        dx_w =  dx_tang2
     else if(surfactant.eq.3)then
        dx_s =  dx_norm
        if (xm(3*i).lt.0.d0 )then
           dx_w = dx_tang 
        else
           dx_w = dx_tang2 
        endif
    else if(surfactant.eq.8)then
        if(plan(i))then
           dx_w =  dx_random
        else
           dx_s =  dx_random
        endif

     endif


!     if(surfactant.eq. 3)then
!        if (xm(3*i).lt.0.d0 )then
!           nrandom = tangential(3*i - 2:3*i)
!        else
!           nrandom = tangential_2(3*i - 2:3*i)
!        endif
!     else  if(surfactant == 2)then
!        if (plan(i))then
!           nrandom = tangential(3*i - 2:3*i)
!        endif
!     endif


!     call RANDOM_NUMBER( nrandom )
!     nrandom = nrandom - normal(3*i - 2:3*i)*dot_product( nrandom, normal(3*i - 2:3*i) )
!     r       = sqrt( dot_product( nrandom, nrandom ) )
!     nrandom = nrandom/r
!    dx = nrandom
!    if(direction.eq.1)then
!         nrandom = tangential(3*i - 2:3*i)
!         dx = normal(3*i - 2:3*i)
!    else if(direction.eq.0)then
!        call RANDOM_NUMBER1( nrandom(1:3) )
!        nrandom = nrandom - 0.5d0
!        r = sqrt( dot_product( nrandom, nrandom ) )
!        nrandom = nrandom/r
!        dx = nrandom
!    endif
     
     call director_tensor( Sbulk, dx_s  , qo_s(:,i) )
     call director_tensor( Sbulk, dx_w  , qo_w(:,i) )

     ! Initial Boundary
     
     q(:,i) = (1.d0 - chi(i) )*qo_w(:,i) + chi(i)*qo_s(:,i)  ! concentration parameter   fs,water=qo_w(:,i) ;  fs,surfactant=qo_s(:,i)
    
  end do
  ! Initial internal
  do i = nn + 1, nodes
     if ( direction == 0 ) then  ! RANDOM
        call RANDOM_NUMBER1( nrandom(1:3) )
        nrandom = nrandom - 0.5d0
        r = sqrt( dot_product( nrandom, nrandom ) )
        nrandom = nrandom/r
        call director_tensor( Sinitial, nrandom, q(:,i) )
       if(wall(i))then
      if (xm(3*i).lt.0.d0 )then
           nrandom = dx_tang
        else
           nrandom = dx_tang2
        endif
 call director_tensor( Sinitial, nrandom, q(:,i) )
 endif

     end if
     if ( direction == 1 ) then  ! UNIFORM
!         if(xm(3*i-2).lt.0.d0)then
          nrandom (1) = 1.d0; nrandom (2) = 1.0d0 ; nrandom(3) = 0.d0
!         else
!          nrandom (1) = -1.d0; nrandom (2) = 0.0d0 ; nrandom(3) = 1.d0
!         endif
 
        r = sqrt( dot_product( nrandom, nrandom ) )
        nrandom = nrandom/r
        call director_tensor( Sinitial, nrandom, q(:,i) )
       
     end if
     if ( direction == 2 ) then   ! RADIAL
        nrandom = xm(3*i - 2:3*i)
        
        r = sqrt( dot_product( nrandom, nrandom ) )
        nrandom = nrandom/r
        
        if(xm(3*i - 2).eq.0.d0.and.xm(3*i - 1).eq.0.d0.and. xm(3*i ).eq.0.d0)then
           nrandom(1)=0.d0 ; nrandom(2)=0.d0 ; nrandom(3)=1.d0
        endif

        call director_tensor( Sinitial, nrandom, q(:,i) )
        
     end if

     if ( direction == 3 ) then   ! BIPOLAR

        nrandom = tangential(3*i - 2:3*i)
        r = sqrt( dot_product( nrandom, nrandom ) )
        nrandom = nrandom/r 

        call director_tensor( Sinitial, nrandom, q(:,i) )
     end if


  end do
 



  ! Checking all the values of the tensor  
!  do i = 1, nodes
!     call eigenval( q(:,i), r, nrandom )
!     r = r*3.d0/2.d0
    
     
!     if ( i <=nn ) SS = Sbulk
!     if ( i > nn ) SS = Sinitial
!     if ( abs(r - SS) > 1.d-10 ) then
!        print *, '--> Problem with initial conditions 1', i, r
!        stop
!     end if
!     r = sqrt( dot_product( nrandom,nrandom ) )
!     if ( abs(r - 1.d0) > 1.d-10 ) then 
!        print *, '--> Problem with initial conditions 2', i, r
!        stop
!     end if
!  end do

  ! Transformation from Q to a

  a = 0.d0
  do i = 1, nodes
     call transformation( q(:,i), a(:,i) )
!     call  check_bound( i, flag )
!     if ( .not. flag ) stop('--> Problem with Initial Conditions 3')
  end do

!  deallocate( q )

end subroutine q_initial
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine director_tensor( s, n, qq )
  use rbf_type
  implicit none
  real(dp), intent(in) :: s
  real(dp), intent(in),  dimension(3) :: n
  real(dp), intent(out), dimension(5) :: qq

  ! this definition uses the SCALAR ORDER PARAMETER
  !   not the maximum eigenvalue 


  qq(1) = s*( n(1)*n(1) - 1.d0/3.d0 )
  qq(2) = s*( n(1)*n(2) )
  qq(3) = s*( n(1)*n(3) )
  qq(4) = s*( n(2)*n(2) - 1.d0/3.d0 )
  qq(5) = s*( n(2)*n(3) )

end subroutine director_tensor
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine eigenval( q, value, vector )
  use rbf_type
  implicit none
  real(dp), intent(in),  dimension(5) :: q
  real(dp), intent(out) :: value
  real(dp), intent(out), dimension(3) :: vector
  integer(i4b) lwork, info,i,j
  real(dp), dimension(3)   :: eigenvalues
  real(dp), dimension(10)  :: work
  real(dp), dimension(3,3) :: a

  lwork  = 10
  ! Matrix form
  a(1,1) = q(1);   a(1,2) = q(2);   a(1,3) = q(3)
  a(2,1) = a(1,2); a(2,2) = q(4);   a(2,3) = q(5)
  a(3,1) = a(1,3); a(3,2) = a(2,3); a(3,3) = - a(1,1) - a(2,2)
  call dsyev( 'V', 'U', 3, a, 3, eigenvalues, work, lwork, info )
  value  = maxval( eigenvalues )
  vector = a(:,3)
 
end subroutine eigenval
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine autovalores_biax( q, value,eta, vector )
  use rbf_type
  implicit none
  real(dp), intent(in),  dimension(5) :: q
  real(dp), intent(out) :: value,eta
  real(dp), intent(out), dimension(3) :: vector
  integer(i4b) lwork, info,i,j
  real(dp), dimension(3)   :: eigenvalues
  real(dp), dimension(10)  :: work
  real(dp), dimension(3,3) :: a

  lwork  = 10
  ! Matrix form
  a(1,1) = q(1);   a(1,2) = q(2);   a(1,3) = q(3)
  a(2,1) = a(1,2); a(2,2) = q(4);   a(2,3) = q(5)
  a(3,1) = a(1,3); a(3,2) = a(2,3); a(3,3) = - a(1,1) - a(2,2)
  call dsyev( 'V', 'U', 3, a, 3, eigenvalues, work, lwork, info )
  value  = maxval( eigenvalues )
  eta = eigenvalues(2)
  vector = a(:,3)
 
end subroutine autovalores_biax
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine transformation( q, a )
  use rbf_type
  implicit none
  real(dp), dimension(5), intent(in) :: q
  real(dp), dimension(5), intent(out) :: a
  
  a = 0.d0
  a(1) = - sqrt(3.d0/2.d0)*( q(1) + q(4) )
  a(2) = sqrt(2.d0)*q(2)
  a(3) = sqrt(2.d0)*q(3)
  a(4) = ( q(1) - q(4) )/sqrt(2.d0)
  a(5) = sqrt(2.d0)*q(5) 

end subroutine transformation
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine inv_transformation( a, q )
  use rbf_type
  implicit none
  real(dp), dimension(5), intent(in) :: a
  real(dp), dimension(5), intent(out) :: q
  real(dp) q33, aux1, aux2

  aux1 = 1.d0/sqrt(6.d0) ; aux2 = 1.d0/sqrt(2.d0)
  q = 0.d0
  q(1) = -a(1)*aux1 + a(4)*aux2
  q(2) =  a(2)*aux2
  q(3) =  a(3)*aux2
  q(4) = -a(1)*aux1 - a(4)*aux2
  q(5) =  a(5)*aux2
  q33 = sqrt(2.d0/3.d0)*a(1)
  q33 = q33 + q(1) + q(4) 
!  if ( abs(q33) > 1.d-10 )stop('Trace is non-zero')

end subroutine inv_transformation
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine initial_da
  use rbf_type
  use variable, only: nn,nip
  use variable, only: da, a
  use variable, only: ix1,dddx,dddy,dddz
  implicit none 
  integer(i4b) i , j, k
  integer(i4b) i1x, i2x, i0, m,ii, i1y, i2y, i1z, i2z
  
  da = 0.d0

  do i = 1,nip
     i0=i+nn
   
     
     i1x=ix1(1,i0)    ! X derivatives
     i2x=ix1(2,i0)

     i1y=ix1(3,i0)    ! Y derivatives
     i2y=ix1(4,i0)

     i1z=ix1(5,i0)    ! Z derivatives
     i2z=ix1(6,i0)
     
!-------------------------------------
      do m=1,5
         
         da(3*i0-2,m)= (a(m,i2x)-a(m,i1x))*dddx
         da(3*i0-1,m)= (a(m,i2y)-a(m,i1y))*dddy
         da(3*i0  ,m)= (a(m,i2z)-a(m,i1z))*dddz
      end do    
  enddo 

         
end subroutine initial_da
!----------------------------------------------------------

!----------------------------------------------------------
!----------------------------------------------------------

subroutine inv_transformation_tensor( a, q )
  use rbf_type
  implicit none
  real(dp), dimension(5), intent(in) :: a
  real(dp), dimension(3,3), intent(out) :: q
  real(dp) q33, aux1, aux2

  aux1 = 1.d0/sqrt(6.d0) ; aux2 = 1.d0/sqrt(2.d0)
  q = 0.d0
  q(1,1) = -a(1)*aux1 + a(4)*aux2
  q(1,2) =  a(2)*aux2
  q(1,3) =  a(3)*aux2
  q(2,2) = -a(1)*aux1 - a(4)*aux2
  q(2,3) =  a(5)*aux2
  q(3,3) =  -q(1,1) - q(2,2)

   
  q(2,1)=q(1,2)
  q(3,1)=q(1,3)
  q(3,2)=q(2,3)

 
end subroutine inv_transformation_tensor
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine transformation_q( q, a )
  use rbf_type
  implicit none
  real(dp), dimension(3,3), intent(in) :: q
  real(dp), dimension(5), intent(out) :: a
  
  a = 0.d0

  a(1) = - sqrt(3.d0/2.d0)*( q(1,1) + q(2,2) )
  a(2) = sqrt(2.d0)*q(1,2)
  a(3) = sqrt(2.d0)*q(1,3)
  a(4) = ( q(1,1) - q(2,2) )/sqrt(2.d0)
  a(5) = sqrt(2.d0)*q(2,3)


end subroutine transformation_q
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
function transition(c)
  use rbf_type
  use variable, only: Lr
  implicit none
  real(dp), intent(in) :: c
  real(dp) transition, k
  k = 100.d0
  transition = 0.5d0*( 1.d0 + tanh(k*( c - Lr ) ))   !!!!!!  transition is the concentration parameter  c=surfactant concentration
end function transition
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine convertion_a_q
  use rbf_type
  use variable, only: a,q,nodes
  implicit none
  integer :: i
  q= 0.d0

  do i=1,nodes 
     call inv_transformation( a(:,i), q(:,i) )
  enddo
end subroutine convertion_a_q
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine convertion_q_a
  use rbf_type
  use variable, only: a,q,nodes
  implicit none
  integer :: i
  a=0.d0
  
  do i=1,nodes 
     call transformation( q(:,i), a(:,i) )
  enddo
end subroutine convertion_q_a
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine q_qtensor(q,q_imc)
  use rbf_type
  implicit none
  real(dp), dimension(3,3), intent(out) :: q_imc
  real(dp), dimension(5), intent(in) :: q

  q_imc(1,1) = q(1)
  q_imc(1,2) = q(2)
  q_imc(1,3) = q(3)
  q_imc(2,2) = q(4)
  q_imc(2,3) = q(5)
  q_imc(3,3) =  -q_imc(1,1) - q_imc(2,2)

   
  q_imc(2,1)=q_imc(1,2)
  q_imc(3,1)=q_imc(1,3)
  q_imc(3,2)=q_imc(2,3)


end subroutine q_qtensor
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine qtensor_q(q_imc,q)
  use rbf_type
  implicit none
  real(dp), dimension(3,3), intent(in) :: q_imc
  real(dp), dimension(5), intent(out) :: q

  q(1) = q_imc(1,1)
  q(2) = q_imc(1,2)
  q(3) = q_imc(1,3)
  q(4) = q_imc(2,2)
  q(5) = q_imc(2,3)
 

end subroutine qtensor_q
!----------------------------------------------------------
!----------------------------------------------------------
!----------------------------------------------------------
subroutine q_qtensor_uni(q,q_imc,i,j)
  use rbf_type
  implicit none
  real(dp) , intent(out) :: q_imc
  real(dp), dimension(5), intent(in) :: q
  integer :: i,j

  if(i.eq.1.and.j.eq.1)  q_imc = q(1) 
  if(i.eq.1.and.j.eq.2)  q_imc = q(2)
  if(i.eq.1.and.j.eq.3)  q_imc = q(3)
  if(i.eq.2.and.j.eq.1)  q_imc = q(2)
  if(i.eq.2.and.j.eq.2)  q_imc = q(4)
  if(i.eq.2.and.j.eq.3)  q_imc = q(5)
  if(i.eq.3.and.j.eq.1)  q_imc = q(3)
  if(i.eq.3.and.j.eq.2)  q_imc = q(5)
  if(i.eq.3.and.j.eq.3)  q_imc = -q(1)-q(4)
 

!  q_imc = q(1)
!  q_imc = q(2)
!  q_imc = q(3)
!  q_imc = q(4)
!  q_imc = q(5)
!  q_imc = -q(1)-q(4)

end subroutine q_qtensor_uni
!----------------------------------------------------------
!----------------------------------------------------------
subroutine local_da(node,index)
  use rbf_type
  use variable, only: nn,dddx,dddy,dddz,wall
  use variable, only: a,da_local
  use variable, only: ix1

  implicit none 
  integer(i4b) i , j, k,node,index
  integer(i4b) i1x, i2x, i0, m,ii, i1y, i2y, i1z, i2z 
  integer ,dimension (7) :: ii_node
  

  ii_node(1)    = node
  ii_node(2:7)  = ix1(1:6,node) 
  da_local = 0.d0  
   m=index

  do i=2,3
     i0 = ii_node(i)

     if (i0.gt.nn) then
        i1x=ix1(1,i0)    ! X derivatives
        i2x=ix1(2,i0)
        da_local(3*i-2,m)= (a(m,i2x)-a(m,i1x))*dddx     
     endif
  enddo

  do i=4,5
     i0 = ii_node(i)
     if (i0.gt.nn) then
        i1y=ix1(3,i0)    ! Y derivatives
        i2y=ix1(4,i0)
        da_local(3*i-1,m)= (a(m,i2y)-a(m,i1y))*dddy
     endif
  enddo

  do i=6,7
     i0 = ii_node(i)
     if (i0.gt.nn) then
        i1z=ix1(5,i0)    ! Z derivatives
        i2z=ix1(6,i0)
        da_local(3*i  ,m)= (a(m,i2z)-a(m,i1z))*dddz
     endif
  enddo

           
end subroutine local_da
